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Abstract A mathematical modeling of Hepatitis C Virus (HCV) dynamics has 
been presented in this paper. The proposed model, which involves four coupled 
ordinary differential equations, describes the interaction of target cells (hepato- 
cytes), infected cells, infectious virions and non-infectious virions. The model takes 
into consideration the addition of ribavirin to interferon therapy and explains the 
dynamics regarding biphasic and triphasic decline of viral load in the model. A 
critical drug efficiency parameter has been defined and it is shown that for efficien- 
cies above this critical value, HCV is eradicated whereas for efficiencies lower this 
critical value, a new steady state for infectious virions is reached, which is lower 
than the previous steady state. 
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1 Introduction 

Hepatitis C is an infectious viral disease caused by Hepatitis C Virus (HCV) and 
its infection route is via the blood. HCV was first identified in 1989, when the 
expression of cDNAs obtained from the blood plasma of a chimpanzee was in- 
duced with hepatitis non-A, non B and was screened with convalescent serum [1]. 
It should be noted that HCV was the first virus to be discovered by molecular 
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biological method and not by previously used virological methods. Accordingly to 
WHO, individuals between 200 to 300 millions worldwide are currently infected 
with HCV [2J. With the discovery of HCV virus, a new dimension on the treatment 
and prevention of liver diseases evolved [3^ . Previously, patients who were thought 
to be suffering from hepatitis non-A or non-B or even alcoholic liver disease, actu- 
ally have hepatitis C. As the dynamics of hepatitis C virus was understood with 
time, it becomes clear that the disease is a major threat for hepatocellular carci- 
noma. 

Since the identification of the hepatitis C virus, great efforts have been made 
to counter-attack it with an antiviral therapy. As a crucial mediator of the innate 
antiviral immune response, interferon-a (IFN-a) was a natural choice for treat- 
ment 4 . With IFN-a therapy, there was a rapid decline in HCV-RNA levels in 
serum but it can achieve only moderate success. With IFN-a monotherapy, that 
is, therapy will IFN-a alone do not achieve clearance of HCV in half of the infected 
individuals suffering from chronic hepatitis C. A major break through came with 
the addition of broad-spectrum antiviral agent ribavirin to IFN-a treatment. With 
combination of pegyleted interferon and ribavirin, a sustained response rates of 
54-56% can be achieved [71[8l[9]. But this combined therapy is expensive and have 
some side effects. Therefore, there is a need for more effective and better tolerated 
therapies for hepatitis C. However, due to difficulty of developing new, potent, 
specific agents against HCV, this combined therapy of pegylated interferon and 
ribavirin is going to remain as the sole therapy in the near future. With more 
better understanding of the mechanism of interferon and ribavirin together with 
improvement in dosing and dose regimen, it is expected that there is going to be 
substantial improvement in the response rates of the patients suffering from hep- 
atitis C. 

Some patients suffering from hepatitis C naturally clear the virus infection 
without any medical intervention but mostly (55-85%) cannot do that and de- 
velop chronic HCV infection [5J. In the standard protocol for the treatment of 
hepatitis C, a patient is given weekly injection of IFN and take ribavirin pills daily 
for the period of treatment 6 . If after six months of therapy, a patient does not 
show any trace of hepatitis C viral load, then the patient is said to have achieved 
sustained virological response (SVR), implying clinically cured. The aim of this 
combined therapy is to reduce the viral load to a minimum, so that a patient 
achieves SVR. If the therapy fails, a patient fails to achieve SVR and suffers from 
chronic hepatitis C. It should be noted that in hepatitis C patients, the rate of 
achievements of SVR depends on the genotype of the particular infecting virus. 
Genotypes f and 4 are usually treated for 48 weeks whereas genotypes 2 and 3 are 
done for 24 weeks with lesser doses of ribavirin lOj. It is observed that 46% of the 
patients with genotype 1 have achieved SVR under this treatment whereas SVR 
rates of 76-82% are observed in patients with HCV genotypes 2 or 3 [11] . 

To understand the dynamics of hepatitis C, specially, to determine the effi- 
ciency of IFN-a as monotherapy and also in combination with ribavirin, mathe- 
matical models have been used extensively, which have provided insights into the 
pathogenesis of HCV in vivo |121I131|T51I16II17] . Other than suggesting predomi- 
nated mechanism(s) of drug action, the models have been successful in explaining 
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the confounding patterns of viral load changes in HCV infected patients undergo- 
ing therapy. 

Neumann et. al [T^], used the basic model of viral dynamics to HCV, which 
assumes a simplified view of HCV infection by taking into account of IFN-a 
monotherapy. The model consists of three populations, namely, uninfected hep- 
atocytes (the target cell) , productively infected hepatocytes and free HCV virions 
and describes the response to interferon therapy. Though the basic model by Neu- 
mann [12 provides insights into the HCV dynamics in vivo but it fails to predict 
the long term response rates observed during combination therapy as it ignores 
the influence of ribavirin. Dixit et. al. ^13^ made an advancement in the basic 
model by taking into account the role of ribavirin action along with interferon. 
The dynamics of model suggests that ribavirin has very little role to play during 
the first phase decline induced by interferon, as ribavirin does not alter the viral 
production. However, ribavirin enhances the second phase slope when interferon 
effectiveness is small, which is in close agreement with experiments [14] . 

One of the drawback of the model by Dixit et. al. [13 is that non-responders 
and patients with triphasic decay patterns are not studied and is restricted to 
biphasic responses only. A modified model, which include the proliferation of un- 
infected and infected cells driven by liver homeostatic mechanisms was studied by 
Dahari et. al. [15.,16] . The model predicts the triphasic decline and succeeds in 
explaining the origins of non-response. In a recent study, Chakrabarty and Joshi 
PJI presented a model where they used deterministic control theory to obtain an 
optimal treatment strategy using interferon and ribavirin. The optimal treatment 
which they obtained succeeded in reducing the levels of viral load and at the same 
time keeping the side effects of the drug to a minimal. Swati DebRoy et al. [18] 
discussed one of the common and life threatening side effect, namely, hemolytic 
anemia, of hepatitis C infection. They used an extension and modification of Neu- 
mann's model [12 to study the effect of combination therapy in the light of anemia 
and succeeded in providing a quantification of the amount of drug a body can tol- 
erate without succumbing to hemolytic anemia. 

In this paper, a mathematical model is proposed depicting the behavior of hep- 
atitis C virus. The idea is to capture the dynamics of the model due to combined 
effect of antiviral drug therapy: interferon and ribavirin. In section 2, formulation 
of the model based of the schematic diagram has been discussed. Qualitative anal- 
ysis of the model has been shown in section 3, which includes determination of 
equilibrium points, positivity and boundedness of the system, local stability anal- 
ysis and persistence of the system. Global stability analysis about the endemic 
equilibrium point has been shown in section 4. Numerical results, capturing the 
dynamics of HCV are discussed in section 5. The paper ends with a discussion. 

2 Mathematical Model 

Adaptive immune responses mediated by T cells are essential in the control of 
HCV and viral clearance. Recent studies in which memory CD4-|- and CD8-I- T 
cells were depleted have confirmed the critical role of these cells in controlling HCV 
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infections [19]. Hence, these cells are the target cells (or hepatocytes). Once in- 
fected with HCV, these target cells becomes infected hepatocytes, which replicates 
the hepatitis C virus. The replicated viruses may be infectious or non-infectious, 
depending on the effectiveness of the drug(s) therapy. Figure 1 gives the schematic 
representation of the above biological scenario. The proposed model is given by 
the following system of coupled ordinary differential equations: 

'^'^ s + rT ( I - '^^] - diT - {1 - cm)aViT (1) 



dt \ k ^ 

^ = {1 - cr,i)aViT - d2l (2) 
''''' -^l-K^)^I-dsVj (3) 



dt \ 2 

dVNI fVr+Vl 



dt 



) (31 - daVNi (4) 



The first equation gives the dynamics of the target cells, that is, hepatocytes. They 
are produced from a source at a constant rate s and at the same time, their growth 
is augmented by a logistic term with an intrinsic growth rate r and carrying capac- 
ity k. The hepatocytes die naturally at a rate di. The logistic term rT (l — ^-^) is 
more realistic in the sense that it incorporates total hepatocytes, both non-infected 
and infected ones. Here 771 denotes the effectiveness (or efficacy) of interferon in 
blocking the release of new virions, crji is fraction of the efficacy (0 < c < 1) and 
hence (1 — crji) gives the ineffectiveness of interferon, which fails to stop the target 
cells getting infected. The number of cells which gets infected is proportional to the 
number of infection virions and available target cells (hepatocytes) with a propor- 
tionality constant a, which explains the term {l — crii)aVjT in (1) and (2). Infected 
cells die at a rate d2. Equations (3) and (4) give the dynamics of virions, namely, 
the infected and non infected ones. The term (-^^^-j^) is the effectiveness of com- 
bined effect of interferon and ribavirin (0<r;r<l,0<r;i <1=>0< < l) 
and hence (l — Hi-illi) is the ineffectiveness of the combined effect, which leads 
to the growth of the infection virions, proportional to the number of infected 
cells (/) with proportionality constant /3 (infection rate). This explain the term 
(1 - aam) pi in (3). The infectious virions die at a rate ds. In (4), the combined 
effect of interferon and ribavirin results in virions (Vni), which are non-infectious 
in nature and also die at a rate ds. System (1-4) has to be analyzed with the 
following initial conditions: T(0) > 0,/(0) > 0, V/(0) > 0,Vni{0) > 0. All system 
parameters are positive. 



3 Qualitative Analysis of the model 
3.1 Positivity and Boundedness 

Theorem 1 The solutions of the system (1-4) are positive for all t > 0. 

Proof Let {T{t), I{t),Vi{t),VNi{t)) be a solution of system (1-4). Consider for 
some time ti > 0,T(ti) = and if ti be the first such time, then < 0. But 

from (1), = s > 0, which is a contradiction. Therefore, T{t) > 0, for all 
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t > 0. 

It is again observed that I{t) > for t > as 

^ > -d2i m = me-"^' (5) 

(6) 

Similarly, (3) and (4) gives 

t 

Vi=(l- ±±211.^ pe-'^'t J I{s)e-''''ds > (7) 



t 

Vj,! = (^^4^) / n^)e-''''ds > 0, (8) 



Therefore, it is concluded that {T{t),I{t),Vi{t),VNi{t)) > for ah t > 0. 
Theorem 2 The solutions of the system (1-4) are bounded. 
Proof Let {T{t), I{t),Vi{t),VNi{t)) be any solution of system. Let 

p{t)=T{t) + I{t) + Vi{t) + VNi{t) (9) 
The derivative of p{t) along the positive solutions of system [1] is given by 
dn rT^ rTl 

-Z=sJ^rT-— — - diT - d2l + PI - dsVi - dsVNi 

<s + rT-^-diT-d2l + pi- daVj - dsVNi 
< s + rT- ^ - Ap{t) 



<-Ap{t) + {s+'^) (10) 

where A = min[di, (^2 — /3), d^] and d2 > p. It follows from (fTOll that 

lim suppm < (s+ ^) = M*. 

Therefore, there exist positive constants M > M* and Ti > such that if t > 
Ti , p{t) < M. This completes the proof. 
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3.2 Equilibria 

System (1-4) has the following non-negative equilibrium points: 



(a) Eo{T, 0,0,0) where T = ^ 



(r-di) + ./(r-di)2 + lp 



, where r > di. We 



assume s < kd\ so that the model is physiologically realistic, that is, T < k. 
(b) E4T*,I*,V;,V^j), where. 



do da 



" (1-C77i)(l- 2421)^/3' 
srR^ + {r-difk{R- 1) 



rR[d2R+{r - di)] 



Vj* 



PI* 



d3 

' Vr+V 



ds 



Here, R, the controlled reproductive number (CRN) [THl is defined as R 
) and the basic reproductive number of the model is calculated as Rq 

Y^. ft may be noted that E^, exists even when _R < 1, provided, srR? + (r 
d\)^k{R — 1) > 0. However, it always exists for _R > 1 



3.3 Local Stability Analysis 

Theorem 3 The disease free equilibrium point Eq is locally asymptotically stable 
if Rq < 1 and unstable if Rq > 1. 



Proof The variational matrix for the equilibrium point Eq is 
f{r-d,)-^ -4 -(l-cr7i)af \ 
-d2 {l-cqi)aT 

(l_Z2n±m)/3 



Vr 







V 



f^)/3 -dj 



The corresponding characteristic equation is 

(^A - (r - di) + (A + d3)(A" -C\ + D)=0 (11) 

where C = -(d2 + ds) and D = [d2d3 - i(l - C77i)(2 - r]r - r]i)a/3f]. For local 
asymptotic stability, all roots of the characteristic equation must be negative or 
have negative real parts. Since, the linear factor gives negative eigenvalues and 

C < 0, for stability D must be positive, that is, < 1. Therefore, the uninfected 
steady state Eq is stable if i?o < 1 and unstable if i?o > 1- 

Theorem 4 The endemic equilibrium point i?* is locally asymptotically stable if 
1<R<3. 
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Proof The variational matrix for the equihbrium point is 

/ srk _ (r-di) (r-di) 



V-2 



k{r-di) JJ 
^Ur-di)fc(d2fl+(r-di))J 

(1- 
( 



R 

-da 



(1 — C77i)afc(r — di) 


\ 


rR 


(1 — CT7i)afc(r — di) 





rR 










-dz) 



The characteristic equation is obtained as 



where 



(Ai + d3)(A3 + AiX" + A2\ + As) 



(r — di) rRs 

da + ds + ^ + 7 — — 

i? (r-di)fc 

(r — di) riis 
R {r-di)k 



(12) 



Aa = (da + da 



((r-di)+dai?) 



(r-di)^ 1 



+ 



rRs 



AxA-2 



As 
As 



dads 
3 



(r-di) 1-^ + 



rRs 



rJ (r-di)fc 



4-1) dadsr+(d^+dl)^^:^ 



+ (da + ds) 
da(r- 



{r-dif 



+ 



2rs 



+ ■ 



r'R's' 



+ 



di 



(r-di) + d2i? 



, , , , , rRs 
+ da(d2 + d3) — 

d2rRs 



i?2 ■ fc ' (r-di)2(fc)2 
(r - di)da + (r - di)ds + 
1 1 



(r-di)^ 



1-1 



+ ■ 



fc((r-di) + d2i?) 



(r-di) (r 
(r - di) + 



■ di) + dai? 
rRs 



+ 



(r-di)fc 



fe(r — di) 



According to Routh-Hurwitz criterion, the necessary and sufficient conditions 
for all the roots of cubic equation (|12p to have negative real parts are Ai > 0, Aa > 
0, As > and Ai Aa - As > 0. Clearly, Ai > 0. Now, Aa > and A3 > if i? > 1 
and AiAa — A3 > provided 1 < < 3. Thus, the endemic equilibrium point 

(T* , /* , Vf* , ) is locally asymptotically stable if the controlled reproductive 
number R satisfies the condition. It should be noted that the endemic equilibrium 
point Et, will be unstable whenever it exists for R < 1. 



3.4 Critical drug efficacy 

In system (1), the effectiveness of interferon and ribavirin are given by the terms 
(1 — crji) and (l — ''''^''^ ). Here, these terms are combined into a single term as 
(1 — C771) (1 — ^ = 1—r^, where 77 represents the overall drug efficiency. Clearly, 

Vi = Vr = 0, that is, drug efficiency is zero before treatment and < 771 < 1,0 < 
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r?r < 1 during antiviral therapy. The stabiUty criteria for uninfected steady state 
{D > 0) is 



d2d3 - i(l - C7?i)(2 -rjr- mhPT > 0, where, T=^ 



(r-di) + -y/(r-di)2 + ^ 



^ (1 - cr]i) [ 1 - ) < - 

1 - < 



2rd2dz 



a/3[(r - + ^(r - di)2fc2 + 4rsfc)] 

This clearly indicates that there exists a point that acts as a point of separation 

between the region of stability for uninfected steady state and the region of stability 
for infected steady state. This point can be termed as a transcritical bifurcation 
point, which is given by 



^_ _ 2rd2d3 

^ ~ aP[{r - di)k + v'(r-di)2fc2+4rsA;)] 

Accordingly, the critical drug efficacy is defined as 

_ ^ 2rd2d3 ^ Tq 

~ aP[{r - di)k + ^(r-di)2fc2+4rsfc)] ~ f 

where Tq = is the number of uninfected hepatocytes in an infected person 

before treatment (obtained by putting 771 = rjr = in T*) and T is the total 
number of hepatocytes in an uninfected individual. 

During antiviral therapy, if 77 > rje, the drug therapy is successful and the 
viral load can be eradicated. On the other hand, if rj < rj^, the viral load and the 
infected cells converge to a new steady state with lower values. 



3.5 Permanence of the system 



Definition 1 A system is said to be permanent if there exist positive constants 
6, A, with 0<5 < A such that 

mini lim infT(t), lim inil{t), lim infVf(t), lim miVNi{t)\>5, 
max< lim supT(f), lim sup7(t), lim snpVi{t), lim anpVNi{t)> < A 

I X— ^ + 00 X— > + oo X— ^+00 X— > + oo I 

for all solutions of the system. 

Theorem 5 The system (1-4) is permanent if (r — di) — (1 — cqijaVi — ^ > 0. 
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Proof Since all the state variables are bounded, it is obvious that I{t) < M** . 
Then (1) can be rewritten as 



dT _^ L T + I 



rT 

diT - (1 - cr]i)aViT <s + rT--— 

k 



(VrT 
V Vk 2 



rk\ rk 

+ T 



rk 



lim supT(t) < { s+ — ] = T (say). 



(3) and (4) gives, 
d{Vi + Vni) 



dt 



Hence, 



Therefore, 



and 



pI-d3{Vi + VNl)<d3 



lim sup(V/ + Vni) < 

t—>-+oo 



lim sup Vj < — - — 

t-)-+oo Cts 



lim sup Vjvj S — r~ 
t->-+oo da 



/3M* 

I3M** 

Vi (say) 
Vni (say). 



(Vr + Vni) • 



Prom (2) we get. 



dl 

c 

So 



— = (1 - cr]i)aViT - d2l < (1 - cqi)ViaT - dil = d2 



(1 - c?7i)Vfar 
d'2 



lim sup/W<il^i;^l^=7. 



Prom the definition of limit superior, again for every e > 0, there exist a time T2 

such that I{t) > (7 - e) for t>T2. 
Again from (3) we get, 



'2-r]i -T]r 



2 - ?7i - ?7r 



) P{I - e) - dsVi 



dz 



{2 -VI- Vrmi - e) 



2d3 



Vi 



Using differential equality we have. 



lim inf Vi(t) > 
t-y+00 



(2 - r?i - rirMl - e) 

2d3 
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Since e > is arbitrarily small, we have 

(2 - 7?l - TJr)fl 



lim inf Vj(t) > 



2d3 



Vi. 



From (4) we get, 

cLVni ^ 
dt 



(r7i+77,)/3(7-e) 



2d3 



Using differential equality we have. 



lim miVNlit) > 

t— J' + OO 

Since e > is arbitrarily small, we have 



Vni 



{m+Vr)l3iI-e) 
2d. 



lim infy^Kt)> ^"^ t7^^' =y.z. 



Again from (1), 



> T 

> T 



T + I 



rl 



diT - {l-cr]i)aViT 
rT' 



{r - di) {I - cr]i)aVi ^ 

{r - di) (1 - crji)aVi — 



{{r-di)-il-cm)aVi-^)--T 

k r 



Therefore, 



if 



lim inf T(t) > 

t— >- + oo 



rl 

(r - di) - (1 - crjijoVi — 



rl 

(r - d\) - (1 - crji)aVi — 



> 0. 



From (2) we get. 



(13) 



^ = {I - crji)aViT - d2l > {1 - c'ni)aViT- d2l = d2[— f!^Llf!^_/]. 

dt — 02 



which implies 



Therefore, from condition p3|) we conclude that T > 0, which ensures that / > 
0, Vi > 0, Vni > 0. This completes the proof of Theorem. 



(1 — crii)aVrT 

lim inilU) > , — - = /. 

d2 
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4 Global Stability of the endemic equilibrium point 



To study the global stability of endemic equilibrium point _E*, the geometric ap- 
proach of Li and Muldowney 20 1 , which guarantees the global stability of endemic 
equilibrium point by obtaining simple sufficient conditions, has been used. 

Theorem 6 Let x f{x) G -R* he a function (class of functions whose deriva- 
tives are continuous) for x in a simply connected domain D C -R^, where 




id f{x) 



/s + rT(l-^) -diT- 
- crji)aViT ■ 



(1 
(1- 



{1 - c7ji)aViT\ 
-d2l 
dsVi 



V 



V 2 , 



131 - d^VN 



J 



Consider the system of differential equations x = f{x) subject to the initial con- 
dition (To, /q, V/o, Vat/o)^ = xo{say) Let x(t, xq) be a solution of the system. The 
system (1-4) has a unique endemic equilibrium point _E* in D and there exits a 
compact absorbing set K G D. It is further assumed that the system (1-4) satisfies 
Bendixson criterion \20^ . that is robust under local perturbations of f at all 
non- equilibrium non-wandering points of the system. Let x — > M(x) fee a 6 x 6 
matrix valued function that is C forx G D. It is also assumed that M~ [x) exists 
and is continuous for x & K . Then the endemic unique equilibrium point is 
globally stable in D if 
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lim sup sup 



XQek 



IJ,{B{x{s,xo)))ds < 



(14) 



where B = MfM + M i^^ M , the value Mf is obtained by replacing each 
entry rUij in M by its directional derivative in the direction of f, \7m*jf and jJ,{B) 
is the Lozinskii measure of B with respect to a vector norm \.\ in i?*, defined by 



KB) 



lim 



\I + hB\ - 1 
h 



(15) 



Proof Since permanence with respect to a set of variables implies persistence with 
respect to the same set of variables, it can be concluded that the system (1-4) is 
persistent. Clearly, the system (1-4) has a unique endemic equilibrium point i?* in 
D and the persistence of the system, together with the boundedness of solutions, 
implies the existence of a compact absorbing set K <Z D [22] . 
The Jacobian matrix J of the system (1-4) is 



I k T 

(1 — cr]i)aVi 
(^ 



V 







rT 
k 

-d2 



-(1 
(1- 



)I3 



- crji )aT 
crii)aT 
-ds 
-ds / 



and the corresponding associated second compound matrix jl^l (for detailed dis- 
cussion of compound matrices, their properties and their relations to differential 
equations, the readers are referred to [23] and [23]) is given by 
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J[2] 



/ ail 







V 



(1 - cr]i)aT (1 - cr]i)aT 

a22 -if 

a22 

(1 — cr]i)aVi — (i2 — ^3 

(1 - c'r]i)aVi 








-(l-cr,i)aT 





-d2 - ds 



then 



(1 



2 /"^ V 2 
J J 





cr)i)aT 

-2d3 



Set the function Q = Q{T, I, Vi, Vni) = diag y\, 1, 1, 1, , 

QfQ-^ = diag (0,0,0,0,1 - - ^) andB = Q/Q-i+Qjl2lQ-i 
Therefore, 



V 



Oil 


(1 


— crji)aT 





(1 - cr]i)aT 








\ 






022 





rT 

k 

















a22 





rTVi 
Ik 


(l-CT)i)aTVj 
J 







(1 


- crji)aVi 





-d2 - dz 















(1 


— crii)al 





ass 


(1 - c?7i)aT 















(- 




aee 


/ 



-Bii -B12 
-B21 -B22 



where 



ail 
a22 
ass 



7^ — d2, 



rT 

k 

rT 

"T 



da, 

i Vj 



Bii = [an] = [- 



rT 
k 



B21 = 0, 0, 0]^ and 



-d2 - d3 + 
aee = -2dz 
^2], B12 = [(1 - crii)aT (1 - cr]i)aT O], 
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Vi' 



B22 



a22 








a22 


cr]i)aVi 





(1 


— cri\)al 









rT 
k 



-d2 - dz 

as 





rTVi 
Ik 







(l-c>?i)aTVr 
J 



(1 - cr]i)aT 
aee 



We now define Lozinskii measure of matrix B as follows: 

< max{gi,g2} (16) 
/x(Bii) + II-B12II and g2 = II-B21II + fi{B22) (||-|| denotes the vector 



where gi 
norm) . 
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It can be shown that (see Appendix) 

^,{B)<j-b (17) 

where 

b=—-f3-{l- crjijaVi - 2(1 - crii)aT - \^ j — 

- max|-d2+ (^^^^^^^) /3,-d3 + (l-c77i)QT} (18) 

for sufficiently large t. Then along each solution {T{t), I{t),Vi{t),VNi{t)) such 
that (T(0),7(0), Vf(0), VAr/(0)) e E we have for t > Z 

Thus boundedness of I{t) and definition of finally gives ^2 < 0. The endemic 
equilibrium point E^,{T* , /*, V/*, Vjv/) is globally stable. 



5 Numerical Results 

System (1-4) has been simulated for viral kinetics during successful antiviral treat- 
ment with the parameter values taken from table 1. Fig. 2 shows the kinetics of 
the viral decline in patients responding to interferon, which is biphasic in nature. 
During the first phase, almost all patients treated with interferon show rapid dose 
dependent viral decline for about 24 to 48 hours |12ll25l[26] . A second phase starts 
after about 48 hours where the viral decline is slow. This dynamics is well captured 
in fig. 2, which shows that the viral load declines to a very low level and eventually 
eradicated, depending on the efficacy of interferon. 

When the effectiveness of interferon is quite large, ribavirin does not have 
significant impact on the second phase decline [27 . This is evident from fig. 3. 
Comparing fig. 2 with fig. 3, it is observed that both of them portray almost iden- 
tical graphs though the efficacy of ribavirin in fig. 2 is zero {ijr = 0) and in fig. 3 is 
0.3 {r]r = 0.3), confirming the fact that ribavirin fails to alter the viral load when 
the efficacy of interferon is effectively large. 

It has been reported that there are few patients who do not respond to ini- 
tial interferon therapy and are termed as null-responders. Re-treatment of nuU- 
responders to a standard interferon regimen shows a considerable first phase de- 
cline, followed by a minor further decline during the second phase and then a 
rebound of HCV is observed [281129]. This dynamics is well portrayed in fig. 4, 
which corresponds to the rebound of the viral load during standard interferon 
regime to the re-treatment of non-responders. 

However, ribavirin enhances the second phase slope if the effectiveness of inter- 
feron is significantly less that 1 (see fig. 5). This interesting behavior is due to the 
fact that interferon initially clears the infectious virions effectively, which results 
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in the first phase decline of the viral load. Ribavirin does not have any effect on 
the first phase decline as it does not alter the viral load. Also, when the effective- 
ness of interferon is large, ribavirin has very little role to play as virion production 
is low but when interferon effectiveness is small, ribavirin renders progeny viri- 
ons non-infectious and enhances the second phase decline. This behavior has been 
confirmed by experimental observations [3011311132] . Fig. 5 is able to capture this 
dynamics with lower interferon efficacy and varying the effectiveness of ribavirin. 
From the figure it is concluded that the patient attains sustained virological re- 
sponses (SVR) during the combined drug therapy, that is, the concentration of 
the viral particles in the blood falls below the detection limit during therapy and 
remains undetected for 24 weeks even after the therapy is stopped. Patients who 
exhibit SVR are generally cured of the HCV infection. 

Several experiments have been conducted to evaluate ribavirin monotherapy 
in the treatment of chronic HCV ^33ii34.i35| . In all these experiments, no virologic 
end-of-treatment-responses (ETR) were observed; either there is a minor decline 
or HCV RNA level remains constant after 3-6 months of ribavirin monotherapy 
compared with pretreatment values [33..34i| . This dynamics is reflected in fig6. 

Fig. 7 shows the viral kinetics after the therapy cessation. A virus resurgence 
post therapy cessation has been simulated for the system (1-4), with different 
drug efficacies (r/i,??^) = ((0.4, 0.6); (0.4, 0.8); (0.4, 0.99)) from time to 14 days 
and then set 771 and rjr to zero for the rest of the simulation. After 14 days, the 
virus resurges to pretreatment levels within (7-7.5) days of post therapy cessation. 
Thus, the model given by system (1-4), predicts resurgences to pretreatment levels 
after cessation of therapy. 

The triphasic decay of viral load for certain parameter values is illustrated in 
fig. 8. As mentioned earlier, viral production and hence viral load decreases due 
to interferon action during the first phase; as a result of which, the production of 
infected cells by new infections falls and the total number of cells decline. Home- 
ostatic mechanisms act to restore the total number of cells by cell proliferation. 
However, since the proliferation term is only with uninfected cells, homeostatic 
mechanisms result predominantly in the proliferation of the uninfected cells. But 
the rate at which the uninfected cells gets infected is much less than the rate at 
which the immune mediated effective killing of virions continues due to interferon 
effectiveness and this results in the second phase of the triphasic decline. As the 
effectiveness of interferon gradually decreases, the role of ribavirin becomes pre- 
dominant and results in further declining of the viral load, which explains the third 
phase of the triphasic decline of the viral load (see fig. 8). 

As there is an increase in the influx rate of new hepatocytes(s), the triphasic 
viral decay disappears and yields a biphasic decay (see fig. 9). It is also to be noted 
that for different drug efficacies close to 1, both the biphasic viral decay (fig. 2) 
and the triphasic viral decay (fig. 8) is close to the death rate (^2) of infected 
hepatocytes (I). 
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6 Conclusion 

Millions of people all over the world are infected with HCV with high levels of 
mortality [36 . The models illustrating the HCV dynamics provide key insights 
into the HCV pathogenesis in vivo and action mechanism of interferon and rib- 
avirin [37 . To study the dynamics of hepatitis C virus, a mathematical model 
comprising of four coupled ordinary differential equations has been proposed to 
understand the role of ribavirin in interferon therapy. The positivity and bounded- 
ness of the system has been established. The basic reproduction number (-Ro) and 
the controlled reproduction number (i?) of the model have been calculated. It is 
observed mathematically that the local stability of the uninfected steady state de- 
pends on basic reproduction number (-Ro) and infected steady state on controlled 
reproduction number {R). The permanence and global stability analysis of the 
model are established. A critical drug efficacy component r]c has been defined and 
characterized. 

Numerical simulation of system (1) shows some interesting results, which mimic 
the dynamics of hepatitis C virus in patients with successful therapy with inter- 
feron and ribavirin. There is a rapid decline in viral load followed by a second 
slower decline (biphasic decline) until the virus becomes undetectable [12.38.139|. 
A triphasic decay [16] of viral load has also been observed in the model with proper 
choice of interferon and ribavirin effectiveness. The critical drug efficacy parameter 
leads to a point (a transcritical bifurcation point) in the model that determines 
the successful antiviral therapy leading to HCV eradication or leading only to a 
partial response. Thus, the HCV dynamics that the model predicts, provide some 
insights on the possible mechanism for the behavior of the viral load observed in 
the clinic [16]. We sincerely hope that the application of this model will give a 
better understanding in the success of anti- viral therapy in HCV infection. 



7 Appendix 

Calculation of Lozinskii measure of matrix B 

Lozinskii measure of matrix B has been defined as follows: 



/i(S) < max{gi,g2} 

where gi = /i(-Bii) + ||Bi2|| and g2 = \\B2i\\ + /i(-B22) ( 
norm) . 



It is easy to compute )i{Bii) 
IIB21I! =/3, then 



do 



gi = - ^ - ^ - d2 + {1 - cr]i)aT 
g2=P + At(S22 = C) 



(20) 

denotes the vector 



I-B12II = (1 — c77i)qT and 



(21) 
(22) 



The matrix [i322]5x5 is partitioned as -622 = C ■ 



Cii C12 
C21 C22 



where Cii = [022] = [-^-f -d3],Ci2 = [O 0],C2i = [0(l-cr?i)aV,00]^ 
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and C22 



022 
-d2 - d3 

(1 — cr)i)al 



rTV, 
Ik 





(l-c?7i)arVf 
I 





(^^^)/3 



055 (1 - cr]i)aT 

066 



Accordingly, 



where 



53 = m(Cii) + I|Ci2| 



rf3 



34 = (1 - cr]i)aVi) + /i(C22) 



(23) 

(24) 
(25) 



To compute /x(C22), the matrix [C22]5x5 is again partitioned as C22 = D 



Dii D12 
D21 -D22 



, where = [022] = [-^-^-ds], D12 = [o _ (i-cm)^TV, ]^ 



Dai = [0 (1 - cr]i)aT 0]^ and 1)22 
Denoting 



-d2 — ds 






155 





(1 - cr]i)aT 

0-66 



fi{D) < max{g5,56} 
where 55 = At(-Dii) + II-D12II and ge = II-D21II + /x(-D22)- 

We have, /i(£'ii) = - ^ - t - t^s, 

||i?i2|| = ^max [r,(i_c,)a] 

= ^[(1 - cm)a] as (1 - cvi)a > f , ||£>2i|| = (1 - crji)aT. Let 



(26) 



D22=E 



E21 E22 



95 



then 



-—---d3 + —[(1 - crii)a\ 



S6 = (1 - c-r]i)aT + /x(D22 = -E) 



(27) 
(28) 



and 



Eli = [-d2-d3],Ei2= [0 0],S21 



V 2 J Vi 



E22 



a.55 (1 - crj])aT 
( '"T''^ )/3 aee 



Again, we define Lozinskii measure of E as follows: 

fi{E) < max{5(7,98} 
where 57 = + Il-E^i2|| and gs = \\E2i\\ +^(£'22). 



(29) 
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We have, /i(^^ii) = -d2 - ds, ||£i2|| = 0, ||i;2i|| = (^^^) ^ and 

^{£22) = max ass + i ^ 1 13, (1 - cr]i)aT + age 

= {7 " W " ^ max + [ ^~\~'^^ ) /3, -d3 + (1 - «?i)aT} 



Therefore 



5-7 = -d2 - da (30) 
58 = < 7 - 7^ - da ^ + 



/ J V 2 y y/ 



+ max |-d2 + ( ^^-^) -^3 + (1 - cr?i)aT| (31) 

So, from equation pO|l . pip and inequahty (|29p . we get 



I Vi ] \ 2 ) Vi 
2 - rjr - r]i' 



+ max|-d2+ ( ^|J-^/3,-d3 + (l-c??i)QT| (32) 

From ((251) and we obtain, 

»^{^4-*} + ^<'--'"^+(^)f 

+ max |-d2 + {^—^^^^) -ds + (1 - cr?i)aTj (33) 



- ' ds (34) 



Also, 



Vi \ 2 J Vi 
I_^ {l-cvi)aTVi (35) 

Then from ((33]) and we get, 

56 < y - ^ + 2(1 - c^,)aT + j — 

+ max |-d2 + + (1 _ c77i)aT| (36) 

From ((27]) and ((35]), we get. 



tT s 

g5 = J + d2 - — - — - ds (37) 
From ((Ml,® and ((23, we get, 

+ max |-d2 + (^~\~'^'^ ) /3. -^3 + (1 - c77i)aT| (38) 
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From ([251) and dSS]), we get, 

54 < y - — + (1 - C77i)aVi + 2(1 - C77i)qT + ^ j — 

+ max |-d2 + ~ ~ ) -ds + (1 - C77i)ar} (39) 
From daSl), (O and ([Ml), we get, 



/i(C^) < y - + (1 - criijaVi + 2(1 - c77i)aT + ^ j — 

+ max |-d2 + { ^~\~'^' ) -ds + (1 - cm)aT^ (40) 
From ((22]) and (gO]), it follows 

32<y-^+/3 + (l- cr;r)ayi + 2(1 - cr?i)aT + (^^^) ^ 

+ max {-da + ( ^~'^^~'" ) -ds + (1 - c77i)aT} (41) 
From (|2ip and psp . we obtain 



/ rT s (1 — cr]i)aVrT , ^ „ ,, ^ 

g^ = --—---\ +(l-C77i)c.T (42) 



Therefore, we get 



where 



m(-B) < y — & for sufficiently large t. (43) 



6=— -/3-(l- c77i)qVi - 2(1 - cT]i)aT - y ^ J — 

- max |-d2 + (^^^^^^) /3. -rfs + (1 - c77i)qt} (44) 
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Table 1 Parameter values used for numerical calculations 



Parameters 


Values 


s (constant rate uninfected hopatocytcs production) 


1.0 cell /ml /day \TB- 


r (proliferation rate) 


1.99/day 14; (OA? - 3.7)/day [40!. 


k (carrying capacity) 


3.6 X W cells/ml 115|. 


di (natural death of hepatocytes) 


0.01/day 


a (rate of infection of uninfected hepatocytes) 


2.25 X W ml /day /virions 15 ; (1.5 - 27.8) X 19^ ' ml / day /virions 1401. 


d2 (natural death rate of infected hepatocytes) 


l.O/day 15 ; (0.43 - 3.1)/day 40 . 


13 (rate at which virions are replicated) 


2.9 virions / cell / day 15 ; (1.2 — 7.9) virions / cell / day |40|. 


^3 (natural death rate of infected and non infected virions) 


6.0/day il5j. 



Note: The effectiveness of the drugs rji and rjr lies between and 1. They are chosen accordingly 
and are mentioned in the figures. The value of c also lies between and 1. 
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Fig. 1 The figure shows the schematic diagram, explaining the dynamics of Hepatitis C virus 
infection. 
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Fig. 2 The figure shows the kinetics of viral decline in patients responding to interferon only, 
which is characterized by biphasic graphs. The effect of ribavirin (rjr) is taken to be zero in 
this case. Three biphasic decline has been shown for various drug efficiencies, of which, two 
of them are greater than the critical drug efficiency rjc. All the parameter values are taken 
from Table 1. 




Fig. 3 The figure shows that ribavirin has no role to play in reducing the viml load, when 
the effectiveness of interferon is large. This is evident as viral load has not been altered when 
compared to figure 2, where effectiveness of ribavirin is zero. All the parameter values are 
taken from Table 1. 
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Fig. 4 The figure sfiows the re-tTeatment of non-responders to a standard interferon regime, 
which shows a first phase decline, followed by further minor decline during the second phase 
and then a rebound of viral load, which are sometimes observed in patients. All the parameter 
values are taken from Table 1. 




Fig. 5 The figure shows the biphasic decline of the viral load when the effectiveness of inter- 
feron is small. All the parameter values are taken from Table 1. 
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Fig. 6 The figure shows the dynamics of ribavirin monotherapy in the treatment of elironic 
HCV. There is very little decline in the viral load, even with high doses of ribavirin. All the 
parameter values are taken from Table 1. 
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Fig. 7 The figure shows the viral kinetics after the therapy cessation. With different drug 
efficiencies, the decline in viral load is observed from to 14 days. After 14 days, the drug 
efficiencies are set to zero and the virus resurges to pretreatment levels within (7-7.5) days of 
post therapy cessation. All the parameter values are taken from Table 1. 
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Fig. 8 T/te figure shows the dynamics of triphasic decline in the viral load. All the parameter 
values are taken from Table 1. 
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Fig. 9 The figure shows the effect of higher influx rates of new hepatocytes, where the triphasic 
decline of viral load changes to biphasic decline. All the parameter values are taken from Table 
1. 



